{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparation of the $\\mu$ opioid receptor with ligand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a complex build system as it has several components, the protein, a sodium ion, the ligand and of course the membrane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from htmd.ui import *\n",
    "from htmd.home import home\n",
    "#get the files\n",
    "shutil.copytree(home()+'/data/mor','/tmp/testmor/pdb')\n",
    "os.chdir('/tmp/testmor')\n",
    "path='./01_prepare/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%ls /tmp/testmor/pdb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#Protein 4dkl is taken from opm\n",
    "\n",
    "topos  = charmm.defaultTopo() + ['pdb/ff.rtf']\n",
    "params = charmm.defaultParam() + ['pdb/ff.prm']\n",
    "prot = Molecule('pdb/4dkl.pdb')\n",
    "prot.filter('protein and noh and chain B or water within 5 of (chain B and protein)')\n",
    "pcenter = np.mean(prot.get('coords','protein'), axis=0)\n",
    "prot = autoSegment(prot, sel='protein') \n",
    "\n",
    "prot = charmm.build(prot, topo=topos, param=params, outdir= path+'prot',ionize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "prot.view()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#Add sodium in the receptor\n",
    "sod = Molecule('pdb/sod.pdb')\n",
    "sod.set('segid','S1')\n",
    "prot.append(sod)\n",
    "\n",
    "#Use a POPC membrane created with vmd and C36\n",
    "memb = Molecule('pdb/membrane80by80C36.pdb')\n",
    "mcenter = np.mean(memb.get('coords'),axis=0)\n",
    "memb.moveBy(pcenter-mcenter)\n",
    "mol = prot.copy()\n",
    "mol.append(memb, collisions=True)  # Append membrane and remove colliding atoms\n",
    "\n",
    "#Add ligand, previously parametrized using gaamp\n",
    "lig = Molecule('pdb/QM-min.pdb') \n",
    "lig.set('segid','L')\n",
    "lcenter = np.mean(lig.get('coords'),axis=0)\n",
    "newlcenter = [np.random.uniform(-10, 10), np.random.uniform(-10, 10),  43]\n",
    "lig.rotateBy(uniformRandomRotation(), lcenter)\n",
    "lig.moveBy(newlcenter - lcenter)\n",
    "mol.append(lig) \n",
    "\n",
    "#Add water\n",
    "coo = mol.get('coords','lipids or protein')\n",
    "m = np.min(coo,axis=0) + [0,0,-5]\n",
    "M = np.max(coo,axis=0) + [0,0,20]\n",
    "mol = solvate(mol, minmax=np.vstack((m,M)))\n",
    "\n",
    "#Build\n",
    "mol = charmm.build(mol, topo=topos, param=params, outdir=os.path.join(path,'build'), saltconc=0.15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "from htmd.protocols.equilibration_v3 import Equilibration\n",
    "from htmd.mdengine.acemd.acemd import GroupRestraint\n",
    "\n",
    "# Use a 10A flat bottom potential to prevent the ligand from diffusing from original position during equilibration\n",
    "width = np.array([10, 10, 10])\n",
    "flatbot = GroupRestraint('segname L and noh', width, [(5, '0ns')])\n",
    "\n",
    "md = Equilibration()\n",
    "md.runtime = 40\n",
    "md.timeunits = 'ns'\n",
    "md.temperature = 300\n",
    "md.restraints = [flatbot] + md.defaultEquilRestraints('20ns')\n",
    "md.useconstantratio = True\n",
    "md.write(os.path.join(path,'build'), os.path.join(path,'equil'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Visualize the flat bottom potential box\n",
    "mol.view('not water')\n",
    "fbcentre = mol.get('coords', sel='segid L').mean(axis=0).squeeze()\n",
    "b = VMDBox(np.vstack((fbcentre - width/2, fbcentre + width/2)).T.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "mdx = LocalGPUQueue()\n",
    "mdx.submit(os.path.join(path, 'equil'))\n",
    "mdx.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the last frame of the equilibration\n",
    "mol = Molecule(os.path.join(path,'equil','structure.psf'))\n",
    "mol.read(os.path.join(path,'equil','output.xtc'))\n",
    "mol.dropFrames(keep=mol.numFrames-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from htmd.protocols.production_v6 import Production\n",
    "from htmd.mdengine.acemd.acemd import GroupRestraint\n",
    "\n",
    "# Apply a flat bottom potential to prevent the ligand from entering from periodic image of the protein\n",
    "width = np.array([70, 70, 60]) \n",
    "\n",
    "# Center the box at residue 218 which is on the upper side of the protein\n",
    "fbcentre = mol.get('coords', sel='protein and resid 218').mean(axis=0).squeeze()\n",
    "flatbot = GroupRestraint('segname L and noh', width, [(5, '0ns')], fbcentre=fbcentre)\n",
    "\n",
    "md = Production()\n",
    "md.runtime = 50\n",
    "md.timeunits = 'ns'\n",
    "md.temperature = 300\n",
    "md.restraints = flatbot\n",
    "md.write(os.path.join(path,'equil'), os.path.join(path,'prod'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "mol.view('not water')\n",
    "b = VMDBox(np.vstack((fbcentre - width/2, fbcentre + width/2)).T.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
